function [allResults, onlyAverageBetas] = extractROImask(subjID, t_thresh)

% PS: If t_thresh is greater than 10, it is intepreted as the desired
% number of voxels that the ROI should have, and the t-value is calculated


allResults = cell(5,5);
allResults{1,1} = subjID;

allResults{1,2} = 'Average t-value';
allResults{1,3} = 'Max t-value';
allResults{1,4} = 'Average Beta';
allResults{1,5} = 'Max Beta';

allResults{2,1} = 'FOBS mask';
allResults{3,1} = 'Adapt mask';
allResults{4,1} = 'Intersection';
allResults{5,1} = 't-threshold';


% Define some basic directories and files
subjectDir = ['/Volumes/cluster/jet/mattar/AdaptID/Analyses/' subjID '/'];
FOBSdir = [subjectDir 'FOBSLoc.feat/'];
FOBSfileZIPED = [FOBSdir 'stats/tstat1.nii.gz'];

% Gunzip FOBS data
FOBSfile_filename = gunzip(FOBSfileZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
FOBSfile_filename = FOBSfile_filename{1};

% Load FOBS data .nii file in MATLAB
FOBSroiStruct = load_untouch_nii(FOBSfile_filename);
FOBSroi = FOBSroiStruct.img;
if t_thresh > 10
    [~,sortIndex] = sort(FOBSroi(:),'descend');
    FOBSroi_mask = zeros(size(FOBSroi));
    FOBSroi_mask(sortIndex(1:t_thresh)) = 1;
else
    FOBSroi_mask = FOBSroi > t_thresh;
end

maskSize = sum(FOBSroi_mask(:));

% Save FFA mask in subject's masks directory
FFAmaskStruct = FOBSroiStruct;
FFAmaskStruct.img = FOBSroi_mask;
if exist([subjectDir 'masks'],'dir') == 0
    mkdir([subjectDir 'masks']);
end
save_untouch_nii(FFAmaskStruct, [subjectDir 'masks/FFAmask.nii']);

% Delete FOBS .nii file
eval(['delete ' FOBSfile_filename]);





% Retrieve the analysis directories
FaceDirectories = struct2cell(dir([subjectDir 'Face*+.feat']));
FaceDirectories = FaceDirectories(1,:)';
numRuns = length(FaceDirectories);

cope_feat = cell(numRuns,1);
varcope_feat = cell(numRuns,1);
tstat_feat = cell(numRuns,1);

FFA_averageT = zeros(numRuns,1);
FFA_maxT = zeros(numRuns,1);
FFA_averageBeta = zeros(numRuns,1);
FFA_maxBeta = zeros(numRuns,1);

for runIndx = 1:numRuns
    % Generate filenames
    copeZIPED = [subjectDir FaceDirectories{runIndx} '/stats/cope3.nii.gz'];
    varcopeZIPED = [subjectDir FaceDirectories{runIndx} '/stats/varcope3.nii.gz'];
    tstatZIPED = [subjectDir FaceDirectories{runIndx} '/stats/tstat3.nii.gz'];
    % Gunzip cope, varcope and tstat
    cope_filename = gunzip(copeZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
    cope_filename = cope_filename{1};
    varcope_filename = gunzip(varcopeZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
    varcope_filename = varcope_filename{1};
    tstat_filename = gunzip(tstatZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
    tstat_filename = tstat_filename{1};
    
    % Load data into matlab
    this_cope = load_untouch_nii(cope_filename);
    this_cope = this_cope.img;
    this_varcope = load_untouch_nii(varcope_filename);
    this_varcope = this_varcope.img;
    this_tstat = load_untouch_nii(tstat_filename);
    this_tstat = this_tstat.img;
    
    % Save maps
    cope_feat{runIndx} = this_cope;
    varcope_feat{runIndx} = this_varcope;
    tstat_feat{runIndx} = this_tstat;
    
    % Delete gunziped files
    eval(['delete ' cope_filename]);
    eval(['delete ' varcope_filename]);
    eval(['delete ' tstat_filename]);
    
    % Mask maps
    maskedBetamap = this_cope .* FOBSroi_mask;
    maskedTmap = this_tstat .* FOBSroi_mask;
    
    % Save important information about this scan
    FFA_averageT(runIndx) = sum(maskedTmap(:))/maskSize;
    FFA_averageBeta(runIndx) = sum(maskedBetamap(:))/maskSize;
    FFA_maxT(runIndx) = max(maskedTmap(:));
    FFA_maxBeta(runIndx) = max(maskedBetamap(:));
    
end

allResults{2,2} = FFA_averageT;
allResults{2,3} = FFA_maxT;
allResults{2,4} = FFA_averageBeta;
allResults{2,5} = FFA_maxBeta;


%% Leave one out approach

Adapt_averageT = zeros(numRuns,1);
Adapt_maxT = zeros(numRuns,1);
Adapt_averageBeta = zeros(numRuns,1);
Adapt_maxBeta = zeros(numRuns,1);

Intersection_averageT = zeros(numRuns,1);
Intersection_maxT = zeros(numRuns,1);
Intersection_averageBeta = zeros(numRuns,1);
Intersection_maxBeta = zeros(numRuns,1);

resulting_t_threshold = zeros(numRuns,1);
numVox_FFAmask_abovethreshold = zeros(numRuns,1);
numVox_Adaptmask_abovethreshold = zeros(numRuns,1);

for leftout = 1:numRuns
    averageover = 1:numRuns;
    averageover = averageover(averageover ~= leftout);
    
    % Generate the t-stat map using the leave-one-out approach
    sumofoneovervarcope = zeros(size(varcope_feat{1}));
    sumofcopeovervarcope = zeros(size(varcope_feat{1}));
    for thisRun = averageover
        sumofoneovervarcope = sumofoneovervarcope + 1./varcope_feat{thisRun};
        sumofcopeovervarcope = sumofcopeovervarcope + cope_feat{thisRun}./varcope_feat{thisRun};
    end
    
    varcope_gfeat = 1./sumofoneovervarcope;
    cope_gfeat = varcope_gfeat .* sumofcopeovervarcope;
    
    tstat_gfeat = cope_gfeat ./ sqrt(varcope_gfeat);
    tstat_gfeat(isnan(tstat_gfeat)) = 0;
    
    
    % Save FaceAdapt mask in subject's masks directory
    FaceAdaptStruct = FOBSroiStruct;
    FaceAdaptStruct.img = tstat_gfeat;
    save_untouch_nii(FaceAdaptStruct, [subjectDir 'masks/FaceAdapt_leaveout' num2str(leftout) '.nii']);
    
    
    % Generate the mask for calculating average effect
    if t_thresh > 10
        [~,sortIndex] = sort(tstat_gfeat(:),'descend');
        FaceAdapt_mask = zeros(size(tstat_gfeat));
        FaceAdapt_mask(sortIndex(1:t_thresh)) = 1;
    else
        FaceAdapt_mask = tstat_gfeat > t_thresh;
    end
    maskSize = sum(FaceAdapt_mask(:));
    
    % Calculate average and max effects within mask
    maskedBetamap = cope_feat{leftout} .* FaceAdapt_mask;
    maskedTmap = tstat_feat{leftout} .* FaceAdapt_mask;
    
    % Save important information about this scan
    Adapt_averageT(leftout) = sum(maskedTmap(:))/maskSize;
    Adapt_averageBeta(leftout) = sum(maskedBetamap(:))/maskSize;
    Adapt_maxT(leftout) = max(maskedTmap(:));
    Adapt_maxBeta(leftout) = max(maskedBetamap(:));
    
    
    % Now, calculate the effects on the intersection of the matrices
    % Generate the mask for calculating average effect
    minFFAMatrix = min(FOBSroi,tstat_gfeat);
    if t_thresh > 10
        [~,sortIndex] = sort(minFFAMatrix(:),'descend');
        intersectionFFA_mask = zeros(size(tstat_gfeat));
        intersectionFFA_mask(sortIndex(1:t_thresh)) = 1;
        resulting_t_threshold(leftout) = minFFAMatrix(sortIndex(t_thresh));
        numVox_FFAmask_abovethreshold(leftout) = sum(FOBSroi(:) > resulting_t_threshold(leftout));
        numVox_Adaptmask_abovethreshold(leftout) = sum(tstat_gfeat(:) > resulting_t_threshold(leftout));
    else
        intersectionFFA_mask = minFFAMatrix > t_thresh;
    end
    maskSize = sum(intersectionFFA_mask(:));
    
    % Calculate average and max effects within mask
    maskedBetamap = cope_feat{leftout} .* intersectionFFA_mask;
    maskedTmap = tstat_feat{leftout} .* intersectionFFA_mask;
    
    % Save important information about this scan
    Intersection_averageT(leftout) = sum(maskedTmap(:))/maskSize;
    Intersection_averageBeta(leftout) = sum(maskedBetamap(:))/maskSize;
    Intersection_maxT(leftout) = max(maskedTmap(:));
    Intersection_maxBeta(leftout) = max(maskedBetamap(:));
    
end

allResults{3,2} = Adapt_averageT;
allResults{3,3} = Adapt_maxT;
allResults{3,4} = Adapt_averageBeta;
allResults{3,5} = Adapt_maxBeta;

allResults{4,2} = Intersection_averageT;
allResults{4,3} = Intersection_maxT;
allResults{4,4} = Intersection_averageBeta;
allResults{4,5} = Intersection_maxBeta;

allResults{5,2} = resulting_t_threshold;
allResults{5,3} = numVox_FFAmask_abovethreshold;
allResults{5,4} = numVox_Adaptmask_abovethreshold;


% Create simplified versions of the general output
onlyAverageBetas = cell(5,(1+numRuns));
onlyAverageBetas{1,1} = subjID;

onlyAverageBetas{2,1} = 'FOBS mask';
onlyAverageBetas{3,1} = 'Adapt mask';
onlyAverageBetas{4,1} = 'Intersection';
onlyAverageBetas{5,1} = 't-threshold';
onlyAverageBetas{6,1} = 'NumVox in FOBSroi';
onlyAverageBetas{7,1} = 'NumVox in ADAPTroi';

for i=1:numRuns
    onlyAverageBetas{1,i+1} = ['Run #' num2str(i)];
    onlyAverageBetas{2,i+1} = allResults{2,4}(i);
    onlyAverageBetas{3,i+1} = allResults{3,4}(i);
    onlyAverageBetas{4,i+1} = allResults{4,4}(i);
    onlyAverageBetas{5,i+1} = allResults{5,2}(i);
    onlyAverageBetas{6,i+1} = allResults{5,3}(i);
    onlyAverageBetas{7,i+1} = allResults{5,4}(i);
end

save(['workspace_' subjID '.mat']);




